Setup
knitr::opts_chunk$set(echo = TRUE, comment = "#>", dpi = 300)
for (f in list.files(here::here("src"), pattern = "R$", full.names = TRUE)) {
source(f)
}
library(rstan)
library(tidybayes)
library(magrittr)
library(tidyverse)
theme_set(theme_classic() + theme(strip.background = element_blank()))
options(mc.cores = 2)
rstan_options(auto_write = TRUE)
1. Linear model: drowning data with Stan
The provided data drowning in the ‘aaltobda’ package contains the number of people who died from drowning each year in Finland 1980–2019. A statistician is going to fit a linear model with Gaussian residual model to these data using time as the predictor and number of drownings as the target variable. She has two objective questions:
- What is the trend of the number of people drowning per year? (We would plot the histogram of the slope of the linear model.)
- What is the prediction for the year 2020? (We would plot the histogram of the posterior predictive distribution for the number of people drowning at \(\tilde{x} = 2020\).)
Corresponding Stan code is provided in Listing 1. However, it is not entirely correct for the problem. First, there are three mistakes. Second, there are no priors defined for the parameters. In Stan, this corresponds to using uniform priors.
a) Find the three mistakes in the code and fix them. Report the original mistakes and your fixes clearly in your report. Include the full corrected Stan code in your report.
- Declaration of
sigmaon line 10 should bereal<lower=0>. - Missing semicolon at the end of line 16.
- On line 19, the prediction on new data does not use the new data in
xpred. This has been changed toreal ypred = normal_rng(alpha + beta*xpred, sigma);.
Below is a copy of the final model. The full Stan file is at models/assignment07-drownings.stan.
data {
int<lower=0> N; // number of data points
vector[N] x; // observation year
vector[N] y; // observation number of drowned
real xpred; // prediction year
}
parameters {
real alpha;
real beta;
real<lower=0> sigma; // fix: 'upper' should be 'lower'
}
transformed parameters {
vector[N] mu = alpha + beta*x;
}
model {
y ~ normal(mu, sigma); // fix: missing semicolor
}
generated quantities {
real ypred = normal_rng(alpha + beta*xpred, sigma); // fix: use `xpred`
}b) Determine a suitable weakly-informative prior \(\text{Normal}(0,\sigma_\beta)\) for the slope \(\beta\). It is very unlikely that the mean number of drownings changes more than 50 % in one year. The approximate historical mean yearly number of drownings is 138. Hence, set \(\sigma_\beta\) so that the following holds for the prior probability for \(\beta\): \(Pr(−69 < \beta < 69) = 0.99\). Determine suitable value for \(\sigma_\beta\) and report the approximate numerical value for it.
#> [1] 0.99225plot_single_hist(x, alpha = 0.5, color = "black") + geom_vline(xintercept = c(-69, 69)) + labs(x = "beta")
c) Using the obtained σβ, add the desired prior in the Stan code.
From some trial and error, it seems that a prior of \(\text{Normal}(0, 26)\) should work. I have added this prior distribution to beta in the model at line 17.
beta ~ normal(0, 26); // prior on `beta`d) In a similar way, add a weakly informative prior for the intercept alpha and explain how you chose the prior.
To use the year directly as the values for \(x\) would lead to a massive value of \(\alpha\) because the values for \(x\) range from 1980 to 2019. Thus, it would be advisable to first center the year, meaning at the prior distribution for \(\alpha\) can be centered around the average of the number of drownings per year and a standard deviation near that of the actual number of drownings.
#> year drownings
#> 1 1980 149
#> 2 1981 127
#> 3 1982 139
#> 4 1983 141
#> 5 1984 122
#> 6 1985 120Therefore, I add the prior \(\text{Normal}(135, 50)\) to \(\alpha\) on line 16.
alpha ~ normal(135, 50); // prior on `alpha`variable_post <- spread_draws(drowning_model, alpha, beta) %>%
pivot_longer(c(alpha, beta), names_to = "variable", values_to = "value")
head(variable_post)
#> # A tibble: 6 × 5
#> .chain .iteration .draw variable value
#> <int> <int> <int> <chr> <dbl>
#> 1 1 1 1 alpha 135.
#> 2 1 1 1 beta -1.12
#> 3 1 2 2 alpha 138.
#> 4 1 2 2 beta -1.57
#> 5 1 3 3 alpha 132.
#> 6 1 3 3 beta -0.691variable_post %>%
ggplot(aes(x = .iteration, y = value, color = factor(.chain))) +
facet_grid(rows = vars(variable), scales = "free_y") +
geom_path(alpha = 0.5) +
scale_x_continuous(expand = expansion(c(0, 0))) +
scale_y_continuous(expand = expansion(c(0.02, 0.02))) +
labs(x = "iteration", y = "value", color = "chain")
variable_post %>%
ggplot(aes(x = value)) +
facet_grid(cols = vars(variable), scales = "free_x") +
geom_histogram(color = "black", alpha = 0.3, bins = 30) +
scale_x_continuous(expand = expansion(c(0.02, 0.02))) +
scale_y_continuous(expand = expansion(c(0, 0.02)))
spread_draws(drowning_model, ypred) %$%
plot_single_hist(ypred, alpha = 0.3, color = "black") +
labs(x = "predicted number of drownings in 2020")
red <- "#C34E51"
bayestestR::describe_posterior(drowning_model, ci = 0.89, test = c()) %>%
as_tibble() %>%
filter(str_detect(Parameter, "mu")) %>%
select(Parameter, Median, CI_low, CI_high) %>%
janitor::clean_names() %>%
mutate(idx = row_number()) %>%
left_join(drowning %>% mutate(idx = row_number()), by = "idx") %>%
ggplot(aes(x = year)) +
geom_point(aes(y = drownings), data = drowning, color = "#4C71B0") +
geom_line(aes(y = median), color = red, size = 1.2) +
geom_smooth(
aes(y = ci_low),
method = "loess",
formula = "y ~ x",
linetype = 2,
se = FALSE,
color = red,
size = 1
) +
geom_smooth(
aes(y = ci_high),
method = "loess",
formula = "y ~ x",
linetype = 2,
se = FALSE,
color = red,
size = 1
) +
labs(x = "year", y = "number of drownings (mean ± 89% CI)")
2. Hierarchical model: factory data with Stan
TODO
#> R version 4.1.1 (2021-08-10)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Big Sur 10.16
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> attached base packages:
#> [1] stats graphics grDevices datasets utils methods
#> [7] base
#>
#> other attached packages:
#> [1] forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7
#> [4] purrr_0.3.4 readr_2.0.1 tidyr_1.1.3
#> [7] tibble_3.1.3 tidyverse_1.3.1 magrittr_2.0.1
#> [10] tidybayes_3.0.1 rstan_2.21.2 ggplot2_3.3.5
#> [13] StanHeaders_2.21.0-7
#>
#> loaded via a namespace (and not attached):
#> [1] nlme_3.1-152 fs_1.5.0 matrixStats_0.61.0
#> [4] lubridate_1.7.10 insight_0.14.4 httr_1.4.2
#> [7] rprojroot_2.0.2 tensorA_0.36.2 tools_4.1.1
#> [10] backports_1.2.1 bslib_0.2.5.1 utf8_1.2.2
#> [13] R6_2.5.0 mgcv_1.8-36 DBI_1.1.1
#> [16] colorspace_2.0-2 ggdist_3.0.0 withr_2.4.2
#> [19] tidyselect_1.1.1 gridExtra_2.3 prettyunits_1.1.1
#> [22] processx_3.5.2 downlit_0.2.1 curl_4.3.2
#> [25] compiler_4.1.1 rvest_1.0.1 cli_3.0.1
#> [28] arrayhelpers_1.1-0 xml2_1.3.2 bayestestR_0.11.0
#> [31] labeling_0.4.2 posterior_1.1.0 sass_0.4.0
#> [34] scales_1.1.1 checkmate_2.0.0 aaltobda_0.3.1
#> [37] callr_3.7.0 digest_0.6.27 rmarkdown_2.10
#> [40] pkgconfig_2.0.3 htmltools_0.5.1.1 highr_0.9
#> [43] dbplyr_2.1.1 rlang_0.4.11 readxl_1.3.1
#> [46] rstudioapi_0.13 jquerylib_0.1.4 farver_2.1.0
#> [49] generics_0.1.0 svUnit_1.0.6 jsonlite_1.7.2
#> [52] distill_1.2 distributional_0.2.2 inline_0.3.19
#> [55] loo_2.4.1 Matrix_1.3-4 Rcpp_1.0.7
#> [58] munsell_0.5.0 fansi_0.5.0 abind_1.4-5
#> [61] lifecycle_1.0.0 stringi_1.7.3 yaml_2.2.1
#> [64] snakecase_0.11.0 pkgbuild_1.2.0 grid_4.1.1
#> [67] parallel_4.1.1 crayon_1.4.1 lattice_0.20-44
#> [70] splines_4.1.1 haven_2.4.3 hms_1.1.0
#> [73] knitr_1.33 ps_1.6.0 pillar_1.6.2
#> [76] codetools_0.2-18 clisymbols_1.2.0 stats4_4.1.1
#> [79] reprex_2.0.1 glue_1.4.2 evaluate_0.14
#> [82] V8_3.4.2 renv_0.14.0 RcppParallel_5.1.4
#> [85] modelr_0.1.8 vctrs_0.3.8 tzdb_0.1.2
#> [88] cellranger_1.1.0 gtable_0.3.0 datawizard_0.2.1
#> [91] assertthat_0.2.1 xfun_0.25 janitor_2.1.0
#> [94] broom_0.7.9 coda_0.19-4 ellipsis_0.3.2
#> [97] here_1.0.1